{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Biome Level statistics for model: Country Level"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import the libraries\n",
    "import ee\n",
    "import pandas as pd\n",
    "import os\n",
    "import numpy as np\n",
    "import random\n",
    "from random import sample\n",
    "import itertools \n",
    "import geopandas as gpd\n",
    "from sklearn.metrics import r2_score\n",
    "from termcolor import colored # this is allocate colour and fonts type for the print title and text\n",
    "from IPython.display import display, HTML"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "#set the working directory of local drive for Grid search result table loading\n",
    "# os.getcwd()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "# initialize the earth engine API\n",
    "ee.Initialize()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1 Load the required composites and images"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load the basic maps that needed for the analysis\n",
    "# load the carbon concentration map\n",
    "carbonConcentration = ee.Image(\"users/leonidmoore/ForestBiomass/Biome_level_Wood_Carbon_Conentration_Map\")\n",
    "# load the root shoot ratio map\n",
    "rootShootRatio = ee.Image(\"users/leonidmoore/ForestBiomass/Root_shoot_ratio_Map\").unmask()\n",
    "# load the two composites tha will be used in the analysis\n",
    "compositeImage =ee.Image(\"users/leonidmoore/ForestBiomass/20200915_Forest_Biomass_Predictors_Image\")\n",
    "compositeImageNew = ee.Image(\"projects/crowtherlab/Composite/CrowtherLab_Composite_30ArcSec\")\n",
    "# load the present forest cover\n",
    "presentForestCover = compositeImage.select('PresentTreeCover')\n",
    "potentialForestCover = ee.Image(\"users/leonidmoore/ForestBiomass/Bastin_et_al_2019_Potential_Forest_Cover_Adjusted\").unmask().rename('PotentialForestCover')\n",
    "\n",
    "# define the present and potential forest cover masks\n",
    "presentMask = presentForestCover.gt(0)\n",
    "potentialMask = potentialForestCover.gt(0)\n",
    "\n",
    "# load the biome layer \n",
    "biomeLayer = compositeImage.select(\"WWF_Biome\")\n",
    "# define a pixel area layer with unit km2\n",
    "pixelAreaMap = ee.Image.pixelArea().divide(10000);\n",
    "stdProj = compositeImage.projection();\n",
    "# define the boundary geography reference\n",
    "unboundedGeo = ee.Geometry.Polygon([-180, 88, 0, 88, 180, 88, 180, -88, 0, -88, -180, -88], None, False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2 Prepare the existing carbon density map: ensemble mean"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load the GS carbon density maps of present biomass\n",
    "carbonDensityGS_Max = ee.Image(\"users/leonidmoore/ForestBiomass/GroundSourcedModel/EnsambledMaps/Predicted_GS1_MaxScaler_Present_density_Ensambled_Mean_20230427\").unmask()\n",
    "carbonDensityGS_Mean = ee.Image(\"users/leonidmoore/ForestBiomass/GroundSourcedModel/EnsambledMaps/Predicted_GS1_MeanScaler_Present_density_Ensambled_Mean_20230427\").unmask()\n",
    "# load the SD carbon density maps of present biomass\n",
    "carbonDensityWK =  ee.Image(\"users/leonidmoore/ForestBiomass/WalkerMap/reprojected_Walker_map_1km\").unmask()\n",
    "# since the SD map from Arnan is biomass density, therefore, we directly transfer thre biomass density into carbon density by multiplying the carbon concentration\n",
    "carbonDensitySD =  ee.Image(\"users/leonidmoore/ForestBiomass/RemoteSensingModel/ESA_CCI_AGB_Map_bias_corrected_1km_2010\").multiply(carbonConcentration).unmask()\n",
    "carbonDensityHM =  ee.Image(\"users/leonidmoore/ForestBiomass/SpawnMap/Spawn_Harmonized_AGB_density_Map_1km\").select('agb').unmask()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3 Prepare the potential carbon density map: ensemble mean"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load the potential carbon density for the GS models\n",
    "potentialDensity_Max1 = ee.Image(\"users/leonidmoore/ForestBiomass/GroundSourcedModel/EnsambledMaps/Predicted_GS1_MaxScaler_Potential_density_Ensambled_Mean_20230427\").unmask()\n",
    "potentialDensity_Max2 = ee.Image(\"users/leonidmoore/ForestBiomass/GroundSourcedModel/EnsambledMaps/Predicted_GS2_MaxScaler_Potential_density_Ensambled_Mean_20230427\").unmask()\n",
    "\n",
    "potentialDensity_Mean1 = ee.Image(\"users/leonidmoore/ForestBiomass/GroundSourcedModel/EnsambledMaps/Predicted_GS1_MeanScaler_Potential_density_Ensambled_Mean_20230427\").unmask()\n",
    "potentialDensity_Mean2 = ee.Image(\"users/leonidmoore/ForestBiomass/GroundSourcedModel/EnsambledMaps/Predicted_GS2_MeanScaler_Potential_density_Ensambled_Mean_20230427\").unmask()\n",
    "\n",
    "# load the potential carbon density maps for SD models\n",
    "potentialDensityWK1 = ee.Image(\"users/nordmannmoore/ForestBiomass/RemoteSensingModel/EnsambleMaps/Predicted_WK1_Potential_density_Ensambled_Mean\").unmask()\n",
    "potentialDensityWK2 = ee.Image(\"users/nordmannmoore/ForestBiomass/RemoteSensingModel/EnsambleMaps/Predicted_WK2_Potential_density_Ensambled_Mean\").unmask()\n",
    "\n",
    "potentialDensitySD1 = ee.Image(\"users/nordmannmoore/ForestBiomass/RemoteSensingModel/EnsambleMaps/Predicted_SD1_Potential_density_Ensambled_Mean\").unmask()\n",
    "potentialDensitySD2 = ee.Image(\"users/nordmannmoore/ForestBiomass/RemoteSensingModel/EnsambleMaps/Predicted_SD2_Potential_density_Ensambled_Mean\").unmask()\n",
    "\n",
    "potentialDensityHM1 = ee.Image(\"users/leonidmoore/ForestBiomass/GroundSourcedModel/EnsambledMaps/Predicted_HM1_Potential_density_Ensambled_Mean\").unmask()\n",
    "potentialDensityHM2 = ee.Image(\"users/nordmannmoore/ForestBiomass/RemoteSensingModel/EnsambleMaps/Predicted_HM2_Potential_density_Ensambled_Mean\").unmask()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define a functiion to do the values adjustment\n",
    "def adjFunc(presentDen,potentialDen):\n",
    "    # check the difference of the two density maps\n",
    "    potentialHigher = potentialDen.subtract(presentDen).gte(0)\n",
    "    potentialLower = potentialDen.subtract(presentDen).lt(0)\n",
    "    # replace the lower potential value by present biomass density value\n",
    "    adjPotentialDen = presentDen.multiply(potentialLower).add(potentialDen.multiply(potentialHigher))\n",
    "    return(adjPotentialDen)\n",
    "\n",
    "# GS max\n",
    "tgbPotentialDensityGS_Max1 = adjFunc(presentDen = carbonDensityGS_Max,potentialDen = potentialDensity_Max1).multiply(rootShootRatio.add(1))\n",
    "tgbPotentialDensityGS_Max2 = adjFunc(presentDen = carbonDensityGS_Max,potentialDen = potentialDensity_Max2).multiply(rootShootRatio.add(1))\n",
    "tgbPresentDensityGS_Max = carbonDensityGS_Max.multiply(rootShootRatio.add(1))\n",
    "# GS Mean \n",
    "tgbPotentialDensityGS_Mean1 = adjFunc(presentDen = carbonDensityGS_Mean,potentialDen = potentialDensity_Mean1).multiply(rootShootRatio.add(1))\n",
    "tgbPotentialDensityGS_Mean2 = adjFunc(presentDen = carbonDensityGS_Mean,potentialDen = potentialDensity_Mean2).multiply(rootShootRatio.add(1))\n",
    "tgbPresentDensityGS_Mean = carbonDensityGS_Mean.multiply(rootShootRatio.add(1))\n",
    "# SD\n",
    "tgbPotentialDensitySD1 = adjFunc(presentDen = carbonDensitySD,potentialDen = potentialDensitySD1).multiply(rootShootRatio.add(1))\n",
    "tgbPotentialDensitySD2 = adjFunc(presentDen = carbonDensitySD,potentialDen = potentialDensitySD2).multiply(rootShootRatio.add(1))\n",
    "tgbPresentDensitySD = carbonDensitySD.multiply(rootShootRatio.add(1))\n",
    "# HM\n",
    "tgbPotentialDensityHM1 = adjFunc(presentDen = carbonDensityHM,potentialDen = potentialDensityHM1).multiply(rootShootRatio.add(1))\n",
    "tgbPotentialDensityHM2 = adjFunc(presentDen = carbonDensityHM,potentialDen = potentialDensityHM2).multiply(rootShootRatio.add(1))\n",
    "tgbPresentDensityHM = carbonDensityHM.multiply(rootShootRatio.add(1))\n",
    "\n",
    "# WK\n",
    "tgbPotentialDensityWK1 = adjFunc(presentDen = carbonDensityWK,potentialDen = potentialDensityWK1).multiply(rootShootRatio.add(1))\n",
    "tgbPotentialDensityWK2 = adjFunc(presentDen = carbonDensityWK,potentialDen = potentialDensityWK2).multiply(rootShootRatio.add(1))\n",
    "tgbPresentDensityWK = carbonDensityWK.multiply(rootShootRatio.add(1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "tgbPresentDensityImage = tgbPresentDensityGS_Max.addBands(tgbPresentDensityGS_Mean).addBands(tgbPresentDensitySD).addBands(tgbPresentDensityHM).addBands(tgbPresentDensityWK)\n",
    "# apply the reducer \n",
    "tgbPresentDensity = tgbPresentDensityImage.reduce(ee.Reducer.mean())\n",
    "\n",
    "# compose the  images into a multiband image\n",
    "tgbPotentialDensityImage = tgbPotentialDensityGS_Max1.addBands(tgbPotentialDensityGS_Max2).addBands(tgbPotentialDensityGS_Mean1).addBands(tgbPotentialDensityGS_Mean2).addBands(tgbPotentialDensitySD1).addBands(tgbPotentialDensitySD2).addBands(tgbPotentialDensityHM1).addBands(tgbPotentialDensityHM2).addBands(tgbPotentialDensityWK1).addBands(tgbPotentialDensityWK2)\n",
    "# apply the reducer \n",
    "tgbPotentialDensity = tgbPotentialDensityImage.reduce(ee.Reducer.mean())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4 Partioning the potential cover into different landuse types"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load all the landuse type layers\n",
    "croplandOrg = ee.Image(\"users/leonidmoore/ForestBiomass/HYDE31/cropland_Percent\").rename('cropland').divide(100).reproject(crs=stdProj);\n",
    "grazingOrg = ee.Image(\"users/leonidmoore/ForestBiomass/HYDE31/grazing_Percent\").rename('grazing').divide(100).reproject(crs=stdProj);\n",
    "pastureOrg = ee.Image(\"users/leonidmoore/ForestBiomass/HYDE31/pasture_Percent\").rename('pasture').divide(100).reproject(crs=stdProj);\n",
    "rangelandOrg = ee.Image(\"users/leonidmoore/ForestBiomass/HYDE31/rangeland_Percent\").rename('rangeland').divide(100).reproject(crs=stdProj);\n",
    "urbanOrg = compositeImage.select(['LandCoverClass_Urban_Builtup']).divide(100).unmask().reproject(crs=stdProj);\n",
    "snowIceOrg = compositeImageNew.select(['ConsensusLandCoverClass_Snow_Ice']).divide(100).unmask().reproject(crs=stdProj);\n",
    "openWaterOrg = compositeImageNew.select(['ConsensusLandCoverClass_Open_Water']).divide(100).unmask().reproject(crs=stdProj);\n",
    "# define the total landcover types\n",
    "sumCover = presentForestCover.add(pastureOrg).add(rangelandOrg).add(croplandOrg).add(urbanOrg).add(openWaterOrg).add(snowIceOrg);\n",
    "oneSubtract = ee.Image(1).subtract(sumCover);\n",
    "freeland = oneSubtract.multiply(oneSubtract.gte(0));\n",
    "# get the scale ratio for pixels with sumCover larger than 1\n",
    "scaleRatio = ee.Image(1).subtract(presentForestCover).divide(sumCover.subtract(presentForestCover)).multiply(oneSubtract.lt(0));\n",
    "# get the ratio of these three disturbed maps\n",
    "pasture = pastureOrg.multiply(scaleRatio).multiply(oneSubtract.lt(0)).add(pastureOrg.multiply(oneSubtract.gte(0))).unmask();\n",
    "rangeland = rangelandOrg.multiply(scaleRatio).multiply(oneSubtract.lt(0)).add(rangelandOrg.multiply(oneSubtract.gte(0))).unmask();\n",
    "cropland = croplandOrg.multiply(scaleRatio).multiply(oneSubtract.lt(0)).add(croplandOrg.multiply(oneSubtract.gte(0))).unmask();\n",
    "urban = urbanOrg.multiply(scaleRatio).multiply(oneSubtract.lt(0)).add(urbanOrg.multiply(oneSubtract.gte(0))).unmask();\n",
    "openWater = openWaterOrg.multiply(scaleRatio).multiply(oneSubtract.lt(0)).add(openWaterOrg.multiply(oneSubtract.gte(0))).unmask();\n",
    "snowIce = snowIceOrg.multiply(scaleRatio).multiply(oneSubtract.lt(0)).add(snowIceOrg.multiply(oneSubtract.gte(0))).unmask();\n",
    "sumTT = presentForestCover.add(pasture).add(rangeland).add(cropland).add(urban).add(freeland).add(openWater).add(snowIce).unmask();\n",
    "\n",
    "effectivePotentialMask = freeland.add(rangeland).add(pasture).add(cropland).add(urban).gt(0);\n",
    "# there are some pixels without any landcover survived but with open water and ice and snow. here we mask these pixels out\n",
    "sumlandCover = pastureOrg.add(rangelandOrg).add(croplandOrg).add(urbanOrg).add(freeland);\n",
    "restorationMap =potentialForestCover.subtract(presentForestCover).mask(effectivePotentialMask).unmask();\n",
    "\n",
    "# sum all these scaled layersv\n",
    "scaledSum = pasture.add(rangeland).add(cropland).add(urban).add(freeland);\n",
    "potentialCoverFinal = restorationMap.add(presentForestCover);\n",
    "# allocate the potential equally to each layer\n",
    "freelandPotentialCover = freeland.divide(scaledSum).multiply(restorationMap).unmask();\n",
    "rangelandPotentialCover = rangeland.divide(scaledSum).multiply(restorationMap).unmask();\n",
    "pasturePotentialCover = pasture.divide(scaledSum).multiply(restorationMap).unmask();\n",
    "croplandPotentialCover = cropland.divide(scaledSum).multiply(restorationMap).unmask();\n",
    "urbanPotentialCover = urban.divide(scaledSum).multiply(restorationMap).unmask();\n",
    "#  allocate the freeland potential in pixels with forest cover larger than 10% to conservation potential\n",
    "freelandForConsevation = freelandPotentialCover.multiply(presentForestCover.gte(0.1)).unmask();\n",
    "maximumPotentialCover = freelandForConsevation.add(presentForestCover);\n",
    "# calucate the reall freeland outside of forest potentialForestCover\n",
    "freelandLeftMap = freelandPotentialCover.subtract(freelandForConsevation).unmask()# the left positive pixels are real freeland pixels"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5 Partioning the biomass potential into different landuse types"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "# calculate the existing carbon, present potential carbon and absolute potential carbon in forests\n",
    "absoluteImage1 = tgbPresentDensity.multiply(pixelAreaMap).multiply(presentMask).divide(1000000000).rename('Present')\n",
    "absoluteImage3 = tgbPotentialDensity.multiply(pixelAreaMap).multiply(potentialCoverFinal.gt(0)).divide(1000000000).rename('AbsolutePotential')\n",
    "# get the sum of the potential covers\n",
    "potentialCoverSum = freelandLeftMap.add(rangelandPotentialCover).add(pasturePotentialCover).add(croplandPotentialCover).add(urbanPotentialCover)\n",
    "\n",
    "trueRestorationPotential = absoluteImage3.subtract(absoluteImage1).multiply(1000000000)\n",
    "ratioPotentialBiomassDensity = absoluteImage1.multiply(potentialCoverFinal.divide(presentForestCover))\n",
    "#  get the real for conservation potential\n",
    "realDensityIncreased = absoluteImage3.subtract(absoluteImage1).mask(absoluteImage3.subtract(ratioPotentialBiomassDensity).gt(0)).unmask()\n",
    "realDensityNotIncreased = absoluteImage3.subtract(absoluteImage1).mask(absoluteImage3.subtract(ratioPotentialBiomassDensity).lte(0)).unmask()\n",
    "trueReforestationPotential = realDensityNotIncreased.add(realDensityIncreased.multiply(ee.Image(1).subtract(presentForestCover.divide(potentialCoverFinal))))\n",
    "\n",
    "conservationPotentialPart1 = realDensityIncreased.multiply(presentForestCover.add(freelandForConsevation).divide(potentialCoverFinal))\n",
    "conservationPotentialPart2 = realDensityNotIncreased.multiply(freelandForConsevation.divide(potentialCoverFinal.subtract(presentForestCover)))\n",
    "\n",
    "absoluteImage2 = conservationPotentialPart1.add(conservationPotentialPart2).add(absoluteImage1).rename('PresentPotential')\n",
    "\n",
    "trueReforestationPotential = absoluteImage3.subtract(absoluteImage2)\n",
    "\n",
    "absoluteImage4 = trueReforestationPotential.multiply(freelandLeftMap.divide(potentialCoverSum)).rename('FreelandPotential')\n",
    "absoluteImage5 = trueReforestationPotential.multiply(rangelandPotentialCover.divide(potentialCoverSum)).rename('RangelandPotential')\n",
    "absoluteImage6 = trueReforestationPotential.multiply(pasturePotentialCover.divide(potentialCoverSum)).rename('PasturePotential')\n",
    "absoluteImage7 = trueReforestationPotential.multiply(croplandPotentialCover.divide(potentialCoverSum)).rename('CroplandPotential')\n",
    "absoluteImage8 = trueReforestationPotential.multiply(urbanPotentialCover.divide(potentialCoverSum)).rename('UrbanPotential')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5.1 Calculate the TGB carbon stock values and write into Google Cloud Storage"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'state': 'READY',\n",
       " 'description': 'Country_Level_Statiscis_of_TGB_export',\n",
       " 'creation_timestamp_ms': 1683721807305,\n",
       " 'update_timestamp_ms': 1683721807305,\n",
       " 'start_timestamp_ms': 0,\n",
       " 'task_type': 'EXPORT_FEATURES',\n",
       " 'id': '7F34PPBJGNBPEBPTCBQ5YNRL',\n",
       " 'name': 'projects/earthengine-legacy/operations/7F34PPBJGNBPEBPTCBQ5YNRL'}"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# read the country borders polygon\n",
    "countryBorders = ee.FeatureCollection(\"users/leonidmoore/ForestBiomass/World_borders_polygon\")#.get('ADM0_NAME')\n",
    "\n",
    "# Stack the absolute biomass layers into an Image.\n",
    "absPotentialImageTGB = absoluteImage1.addBands(absoluteImage2).addBands(absoluteImage3).addBands(absoluteImage4).addBands(absoluteImage5).addBands(absoluteImage6).addBands(absoluteImage7).addBands(absoluteImage8)\n",
    "# define the function to do the country level statistics which could be applied by map      \n",
    "def countryLevelStat(fc):\n",
    "    maskedImg = absPotentialImageTGB.clip(fc.geometry())\n",
    "    output = maskedImg.reduceRegion(reducer= ee.Reducer.sum(),\n",
    "                                    geometry= unboundedGeo,\n",
    "                                    crs='EPSG:4326',\n",
    "                                    crsTransform=[0.008333333333333333,0,-180,0,-0.008333333333333333,90],\n",
    "                                    maxPixels= 1e13)\n",
    "    # generat the eedictionary contains the counrty or region names\n",
    "    countryRegionName = ee.Feature(fc).get('NAME')\n",
    "    fcOut = ee.Feature(ee.Geometry.Point([0,0])).set(output).set('Name',countryRegionName)\n",
    "    return fcOut\n",
    "\n",
    "\n",
    "\n",
    "statisticTable = countryBorders.map(algorithm = countryLevelStat)\n",
    "countryStat = ee.batch.Export.table.toCloudStorage(collection = statisticTable, \n",
    "                                                   description = 'Country_Level_Statiscis_of_TGB_export',\n",
    "                                                   bucket = \"crowtherlab_gcsb_lidong\", #replace it by your own bucket\n",
    "                                                   fileNamePrefix = 'ForestBiomassExport/Country_Level_Statistics_of_TGB_20230427',\n",
    "                                                   fileFormat = \"CSV\")\n",
    "\n",
    "\n",
    "# start the export task\n",
    "countryStat.start()\n",
    "# show the task status\n",
    "countryStat.status()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5.2 Calculate the PGB carbon stock values and write into Google Cloud Storage"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'state': 'READY',\n",
       " 'description': 'Country_Level_Statiscis_of_PGB_export',\n",
       " 'creation_timestamp_ms': 1683721808537,\n",
       " 'update_timestamp_ms': 1683721808537,\n",
       " 'start_timestamp_ms': 0,\n",
       " 'task_type': 'EXPORT_FEATURES',\n",
       " 'id': 'PEJDVXVJCYVMXHAYCGGARVSS',\n",
       " 'name': 'projects/earthengine-legacy/operations/PEJDVXVJCYVMXHAYCGGARVSS'}"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "deadWoodLitterRatio = ee.Image(\"users/leonidmoore/ForestBiomass/DeadWoodLitter/DeadWood_Litter_Ratio_Map\").unmask()\n",
    "\n",
    "# multiply the dead wood and ratio to get the PGB for each country\n",
    "absPotentialImagePGB = absPotentialImageTGB.multiply(deadWoodLitterRatio)\n",
    "# define the function to do the biome level statistics which could be applied by map      \n",
    "def countryLevelStat(fc):\n",
    "    maskedImg = absPotentialImagePGB.clip(fc.geometry())\n",
    "    output = maskedImg.reduceRegion(reducer= ee.Reducer.sum(),\n",
    "                                    geometry= unboundedGeo,\n",
    "                                    crs='EPSG:4326',\n",
    "                                    crsTransform=[0.008333333333333333,0,-180,0,-0.008333333333333333,90],\n",
    "                                    maxPixels= 1e13)\n",
    "    # generat the eedictionary contains the counrty or region names\n",
    "    countryRegionName = ee.Feature(fc).get('NAME')\n",
    "    fcOut = ee.Feature(ee.Geometry.Point([0,0])).set(output).set('Name',countryRegionName)\n",
    "    return fcOut\n",
    "\n",
    "\n",
    "statisticTable = countryBorders.map(algorithm = countryLevelStat)\n",
    "countryPGB = ee.batch.Export.table.toCloudStorage(collection = statisticTable, \n",
    "                                                  description = 'Country_Level_Statiscis_of_PGB_export',\n",
    "                                                  bucket = \"crowtherlab_gcsb_lidong\", #replace it by your own bucket\n",
    "                                                  fileNamePrefix = 'ForestBiomassExport/Country_Level_Statistics_of_PGB_20230427',\n",
    "                                                  fileFormat = \"CSV\")\n",
    "\n",
    "\n",
    "# start the export task\n",
    "countryPGB.start()\n",
    "# show the task status\n",
    "countryPGB.status()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5.3 Calculate the PGB carbon stock values and write into Google Cloud Storage"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load the carbon density layers\n",
    "SandermannCarbonDiff = ee.Image(\"users/leonidmoore/ForestBiomass/SoilOrganicCarbonModel/SOCS_0_200cm_Diff_1km_Present_subtract_NoLU\").unmask()\n",
    "SandermannCarbonPresent = ee.Image(\"users/leonidmoore/ForestBiomass/SoilOrganicCarbonModel/SOCS_0_200cm_1km_Present\").unmask()\n",
    "\n",
    "# mask the diffrence layer\n",
    "SandermannCarbonLoss = SandermannCarbonDiff.multiply(SandermannCarbonDiff.gt(0)).multiply(potentialForestCover.gte(0.1))\n",
    "\n",
    "# calculate the existing carbon, present potential carbon and absolute potential carbon in forests\n",
    "# absoluteImage1 = SandermannCarbonLoss.multiply(pixelAreaMap).multiply(presentForestCover).divide(1000000000).mask(potentialMask)\n",
    "absoluteImage2 = SandermannCarbonLoss.multiply(pixelAreaMap).multiply(maximumPotentialCover).divide(1000000000).mask(potentialMask).rename('ConservationPotential')\n",
    "absoluteImage3 = SandermannCarbonLoss.multiply(pixelAreaMap).multiply(potentialCoverFinal).divide(1000000000).mask(potentialMask).rename('AbsolutePotential')\n",
    "absoluteImage4 = SandermannCarbonLoss.multiply(pixelAreaMap).multiply(freelandLeftMap).divide(1000000000).mask(potentialMask).rename('FreelandPotential')\n",
    "absoluteImage5 = SandermannCarbonLoss.multiply(pixelAreaMap).multiply(rangelandPotentialCover).divide(1000000000).mask(potentialMask).rename('RangelandPotential')\n",
    "absoluteImage6 = SandermannCarbonLoss.multiply(pixelAreaMap).multiply(pasturePotentialCover).divide(1000000000).mask(potentialMask).rename('PasturePotential')\n",
    "absoluteImage7 = SandermannCarbonLoss.multiply(pixelAreaMap).multiply(croplandPotentialCover).divide(1000000000).mask(potentialMask).rename('CroplandPotential')\n",
    "absoluteImage8 = SandermannCarbonLoss.multiply(pixelAreaMap).multiply(urbanPotentialCover).divide(1000000000).mask(potentialMask).rename('UrbanPotential')\n",
    "# absoluteImage9 = SandermannCarbonLoss.multiply(pixelAreaMap).multiply(freelandForConsevation).divide(1000000000).mask(potentialMask)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'state': 'READY',\n",
       " 'description': 'Country_Level_Statiscis_SoilCarbon_export',\n",
       " 'creation_timestamp_ms': 1683721810018,\n",
       " 'update_timestamp_ms': 1683721810018,\n",
       " 'start_timestamp_ms': 0,\n",
       " 'task_type': 'EXPORT_FEATURES',\n",
       " 'id': 'ZIBVHQH4WOW5MTTHESM7Z3Y2',\n",
       " 'name': 'projects/earthengine-legacy/operations/ZIBVHQH4WOW5MTTHESM7Z3Y2'}"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Stack the absolute biomass layers into an Image.\n",
    "absPotentialImageSoil = absoluteImage2.addBands(absoluteImage3).addBands(absoluteImage4).addBands(absoluteImage5).addBands(absoluteImage6).addBands(absoluteImage7).addBands(absoluteImage8)\n",
    "# read the country borders polygon\n",
    "countryBorders = ee.FeatureCollection(\"users/leonidmoore/ForestBiomass/World_borders_polygon\")#.get('ADM0_NAME')\n",
    "\n",
    "# define the function to do the country level statistics which could be applied by map      \n",
    "def countryLevelStat(fc):\n",
    "    maskedImg = absPotentialImageSoil.clip(fc.geometry())\n",
    "    output = maskedImg.reduceRegion(reducer= ee.Reducer.sum(),\n",
    "                                    geometry= unboundedGeo,\n",
    "                                    crs='EPSG:4326',\n",
    "                                    crsTransform=[0.008333333333333333,0,-180,0,-0.008333333333333333,90],\n",
    "                                    maxPixels= 1e13)\n",
    "    # generat the eedictionary contains the counrty or region names\n",
    "    countryRegionName = ee.Feature(fc).get('NAME')\n",
    "    fcOut = ee.Feature(ee.Geometry.Point([0,0])).set(output).set('Name',countryRegionName)\n",
    "    return fcOut\n",
    "\n",
    "\n",
    "\n",
    "statisticTable = countryBorders.map(algorithm = countryLevelStat)\n",
    "countryStat = ee.batch.Export.table.toCloudStorage(collection = statisticTable, \n",
    "                                                   description = 'Country_Level_Statiscis_SoilCarbon_export',\n",
    "                                                   bucket = \"crowtherlab_gcsb_lidong\", #replace it by your own bucket\n",
    "                                                   fileNamePrefix = 'ForestBiomassExport/Country_Level_Statistics_of_SoilCarbon_20230427',\n",
    "                                                   fileFormat = \"CSV\")\n",
    "\n",
    "\n",
    "# start the export task\n",
    "countryStat.start()\n",
    "# show the task status\n",
    "countryStat.status()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# All those tables were donloaded from cloud storage to local folder 'Data/BiomeLevelStatistics/StatisticsForModels'"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
